Complex-tensor theory of simple smectics

Matter self-assembling into layers generates unique properties, including structures of stacked surfaces, directed transport, and compact area maximization that can be highly functionalized in biology and technology. Smectics represent the paradigm of such lamellar materials — they are a state between fluids and solids, characterized by both orientational and partial positional ordering in one layering direction, making them notoriously difficult to model, particularly in confining geometries. We propose a complex tensor order parameter to describe the local degree of lamellar ordering, layer displacement and orientation of the layers for simple, lamellar smectics. The theory accounts for both dislocations and disclinations, by regularizing singularities within defect cores and so remaining continuous everywhere. The ability to describe disclinations and dislocation allows this theory to simulate arrested configurations and inclusion-induced local ordering. This tensorial theory for simple smectics considerably simplifies numerics, facilitating studies on the mesoscopic structure of topologically complex systems.

Layered materials are key components in many technological, biological and fluidic systems. Graphene, MXene and other two-dimensional materials are composed of atomistically thin sheets with only weak out-of-plane bonding 1 . Moreover, layering is a widespread strategy for increasing the surface area of organs, as in the friction ridges that make up fingerprints (Fig. 1a), convolutions within the cortex of human brains 2 , and the enormous surface area of intestinal villi 3 . At the subcellular level, the Golgi apparatus, rough endoplasmic reticulum, and crista in mitochondria all possess many folds and creases. Layered fluids include liquid crystalline cholesterics which form "pseudolayers" via their helicity 4 and, most plainly, smectics composed of stacks of orientationally aligned molecules that maintain fluidic disorder within each layer 5 .
However, the very properties that make these lamellar phases so interesting also contrive to make them challenging to model. Any order parameter that contends to describe simple smectics, i.e., lamellar phases without underlying nematic order, must explicitly include three pieces of information: (i) the degree of ordering, (ii) dilation/contraction of layers' spacing and (iii) the direction/bending of layers, either through its instantaneous value or gradients. Since the 1970s, it has thus been recognized that a complex scalar order parameter can be employed to incorporate this information 20 , which leads to a strong analogy to superfluids and superconductors [21][22][23][24] . However, such a complex order parameter faces many difficulties [25][26][27] . Fundamentally, these issues arise because complex scalar order parameters assume a global layer normal direction is known and that the layers only vary about this well-defined direction.
To circumvent this shortcoming, theoreticians have modeled smectic-A liquid crystals by coupling models of smectic ordering to the nematic order parameters for the orientation of rod-like molecules [28][29][30][31] . While sensible for smectic liquid crystals with inner-layer nematic ordering, this does not fundamentally solve the issue for many non-nematic lamellar materials, which have their own intrinsic rotational symmetry breaking. Furthermore, these approaches model smectics at the microscopic level, resolving each layer at length scales comparable to the molecular size. This makes such models computationally expensive for mesoscale and hydrodynamic modeling.
Recent studies of defect annihilation in block copolymer films 32,33 and confined smectic colloidal liquid crystals showing unique pairs of quarter-charge topological defects connected by domain wall bridges 34,35 motivate the need for alternative theoretical descriptions for simple smectics that allow simulations to tackle more topologically complex structures without relying on microscale models. To create a more general theory for layered materials, we consider "simple" smectics, an idealization of purely lamellar materials that do not simultaneously require a strong degree of orientational alignment of mesogens. We seek a mesoscopic order parameter for such simple smectics that can incorporate both the smectic ordering and layer orientation, and describe our theory as 'mesoscopic' because it neglects details of the material structure on molecular (microscopic) length scales, but is able to describe defects as well as the bulk (macroscale) ordering. Here, we propose a formalism to model defect-rich lamellar systems through a complex tensor order parameter that describes the local degree of lamellar ordering, layer displacement and layer normal orientation.
A lamellar system is a periodic variation of the density in a single spontaneously symmetry-broken direction. Traditionally, smectic ordering is approximated by expanding the molecular density at each point r, as ρ r, t ð Þ= P 1 m = À1 ψ m e imqÁr ≈ ρ 0 + 2Re Ψ ½ , where ρ 0 is the mean and the microscopic density variation is defined to be the real part of a plane wave with Ψ r, t ð Þ= |ψ|e i qÁr + φ 0 ð Þ . Commonly employed in Landau free energy expansions 36,37 , the modulus |ψ r, t ð Þ| quantifies the extent of layering, φ 0 is the phase at the origin, and q r, t ð Þ is the wave vector. In the absence of any deformations, the wave vector q 0 = q 0 N points in the direction of the layer normal N r, t ð Þand describes equally spaced layers of thickness 2π/q 0 . If the layers are locally displaced from their ground state by a layer displacement field u r, t ð Þthat varies slowly on large length and time scales, the argument of the plane wave can be written as Φ r, t ð Þ= q r, t ð ÞÁr≈q 0 Á r + u ð Þ. However, the layer displacement must also be parallel to N by definition since the simple smectic is isotropic within layers. Thus, the phase ϕ r, t ð Þ q 0 Á u encodes the extent of layer displacement in units of the unperturbed wave number q 0 . Variation of ϕ r, t ð Þ indicates lamellar compression/dilation deformations (hereafter referred to jointly as compression). The symmetries of simple smectics can be described as N → − N, Φ → Φ + m2π for m 2 Z, and Φ → − Φ. This last equivalence relation is identical to Ψ → Ψ * . We note that equations such as the Swift-Hohenberg 38 , phase field crystals [39][40][41][42][43] and density functional theory [44][45][46][47] can model striped phases via a continuous field approach that goes beyond the traditional simplistic sinusoidal approximation. However, these all must resolve individual layers, necessarily limiting their applicability to short times and length scales. Though traditionally the microscopic variation is approximated as the lowest Fourier component of the density, the imperative conclusion is that the (i) degree of layering, (ii) layer displacement, and (iii) layer normal direction are the crucial aspects of the lamellar ordering in the continuum limit that vary on mesoscopic scales, regardless of the microscopic details.
Though elegant and economical, the traditional formalism has known shortcomings 26 . These can be stated in a number of equivalent ways: (i) Only the variation of the density ∼ Re Ψ ½ = |ψ| cosðΦÞ is physical, and so the wave function is not truly a single-valued function of position 23,48 but rather both Ψ and Ψ * appear simultaneously 27 . (ii) This can be stated as Φ is not an element of the unit circle S 1 but rather of the orbifold S 1 =Z 2 25,49,50 . (iii) The plane wave is a linear function of the layer normal N, which does not faithfully reflect the apolar symmetry of the layering. While the phase ϕ r, t ð Þ= ðq 0 NÞ Á ðuNÞ is an even function of the layer normal, the first term of Φ r, t ð Þ= q 0 N Á r + ϕ is not. Commonly, the layer normal is defined as a vector via the gradients a Single defects marked on a high resolution photo of a fingerprint. +1/2 disclination, −1/2 disclination and edge dislocation marked by red cross, yellow trilateral and red circle, respectively. b Schematics of (left) +1/2 disclination, (center) −1/2 disclination and (right) dislocations. c-k Simulations of defects for model parameters A = − 1 (lamellar state), C = 2 and κ 2 = 0.75 in circular domains with boundary conditions requiring single defects. Columns present three defect types: (first column) + 1/2 disclination; (second) − 1/2 disclination; (third) edge dislocation. Rows show plots of: c-e layer visualization Re Ψ ½ ; f-h modulus |ψ| with N overlaid; i-k phase ϕ with N.
∇Φ=|∇Φ|, but this is only self-consistent when variations of ϕ are negligible. By explicitly considering microscopic density variation on the scale of individual layers 26,31,41 , these issues can be avoided with computational cost, but a hydrodynamic-scale theory that both circumvents these issues and considers only mesoscopic variations of the lamellar order has not previously been established.
In contrast, hydrodynamic-scale nematic theory uses the tensor Q = S n n À δ=d À Á to collect both the scalar order parameter S and apolar director n into a single order parameter, for dimensionality d and identity matrix δ 51 . The nematic order parameter simultaneously describes the extent of phase ordering and local direction of broken symmetry in an arbitrary reference frame. Thus, both the bulk and deformation free energy densities can written in terms of Q r, t ð Þ. Practically, the introduction of Q has enabled numerical simulations of confined nematics 52,53 , colloidal liquid crystals [54][55][56] and active fluids 57-59 , by treating defects as locally disordered cores, rather than singularities. Explicitly, the tensor Q regularizes these singularities in n and remains continuous at the center of defects. In this manuscript, we show how a similar tensorial order parameter can be constructed for a hydrodynamic-scale description of simple smectics that naturally respects the apolar symmetry of layering, the extent of layering and the relative layer displacement.

Results
In this article, we propose a tensorial order parameter field for simple, lamellar smectics. The tensor E r, t ð Þ is complex, symmetric, traceless and globally gauge invariant. It incorporates the extent of layering and relative layer displacement, as well as the layer normal orientation. It encompasses the advantages Q-tensor formalism provides to nematics but for simple smectics. Here, we exclusively consider the simplest smectics, with only lamellar broken translational symmetry and layernormal broken rotational symmetry. Such simple lamellar ordering is observed in systems with little-to-no nematic ordering, such as the striped phases of short-range attraction and long-range repulsion colloidal systems 44,60,61 , ionic liquid crystals (which commonly transition directly from isotropic to smectic phases) 62 and the flat plate ("lasagna") phase of nuclear pasta 63 . We focus solely on this idealistic premise for lamellae to demonstrate the fundamental suitability of E for avoiding the phase ambiguity and to show its utility in simulating confining geometries.
We consider only the mesoscopic aspects of simple smectics by rearranging the wave function as Ψ = |ψ|e iϕ e iq 0 Ár . In this form, the wave function possesses two distinct factors: (i) The local microscopic ground state wave e iq 0 Ár and (ii) the mesoscopic complex order parameter ψ r, t ð Þ= |ψ|e iϕ . Both the extent of layering |ψ| r, t ð Þ= ψψ * À Á 1=2 and the phase ϕ r, t ð Þ are assumed to vary slowly over mesoscopic times and length scales. Being even in the apolar layer normal N r, t ð Þ, the phase does not possess the same ambiguity as Φ. A change of ϕ by 2π indicates that the layer displacement u has increased by a full layer thickness, which necessitates the existence of an additional layer in this mesoscopic model.
To account for the apolar symmetry of the layer normal, the smectic tensorial order parameter E r, t ð Þ must contain the dyadic square of N, making E symmetric. Furthermore, the absence of preferential directions within planar layers indicates local rotations about N are arbitrary. A traceless order parameter ensures linear terms do not contribute to the bulk free energy. Based on these considerations, we propose the uniaxial complex-tensorial smectic order parameter The scalar order parameter ψ r, t ð Þ= |ψ|e iϕ 2 C is the largest eigenvalue of E (with |ψ| 2 R and ϕ ∈ S 1 ) and the layer normal N r, t ð Þ 2 R d is the associated eigenvector. The tensor describes mesoscopic variations of lamellar order in the hydrodynamic limit, treating lengths on the scale of individual layers as microscopic. It is symmetric, traceless and globally gauge invariant (under E → e iθ E for arbitrary θ); furthermore, it allows N → − N and, by using ϕ, avoids interpretation issues stemming from the double-valued nature of Φ. This limits its consideration only to mesoscopic variations in layer displacement. While a generic second-rank tensor has d 2 elements, these constraints ensure E only possesses d + 1 degrees of freedom, representing the extent of layering, layer displacement, and unit layer normal. As in nematic Qtheory, local biaxiality is possible in E in 3D but biaxiality at equilibrium requires higher order terms in the bulk free energy 64 (see Methods). In describing these slowly varying aspects of a lamellar system without reference to any specific form of the microscopic density variation, E acts as a mesoscopic order parameter that is not based on or limited to any particular molecular details or assumptions about the layer structure.
Though smectics and other lamellae have been modeled from many perspectives 12,33,65,66 , we consider a Landau free energy expansion. The total free energy density f is the sum of bulk and two deformation (compression and curvature) terms. All contributions to the free energy must be real, requiring pairings of E and its complex conjugate E * . Furthermore, the free energy should not depend on E in a manner that is equivalent to a direct dependence on the phase, since it can be globally shifted.

Bulk
Since E is traceless, the bulk smectic free energy density can be written where C > 0, and Einstein summation convention is adopted. Lamellar order is established when A < 0, but the fluid is isotropic when A > 0. The bulk free energy does not depend on phase or layer normal, but only on |ψ|. This form is consistent with scalar-based bulk free energies 36,37,67 (see Methods). In the mean-field limit, Eq. (2) predicts a second order phase transition.

Compression
Lamellae possess two deformation modes: (i) compression, and (ii) curvature of the layers. We consider first compression free energies, which involve derivatives of the tensor order parameter. The simplest such term is E ij,k E * ij,k , where k denotes the direction of the gradient. Additional real terms could be constructed through combinations of similar forms, which would allow different deformation modes to possess differing elastic modulii. For clarity, we make a one-constant approximation where b 1 is a layer compression elastic constant. Equation (3) accommodates first order distortions of the layer normal, as well as contributions due to gradients of the complex amplitude ψ. In a vectorbased model, defects would create topological singularities in the layer normal field, making gradient terms in Eq. (3) irregular; however, E regularizes the singularities and Eq. (3) is continuous.

Curvature
Distortions from uniformly aligned layers come with a free energy cost, akin to a membrane curvature free energy density. We again make a one-constant approximation and keep only the simplest term where b 2 is a bending modulus.
Article https://doi.org/10.1038/s41467-023-36506-z By inserting the eigenvalue and associated eigenvector via Eq. (1) into the free energy contributions (Eqs. (2)-(4)), one can directly compare E-theory to existing models of smectics (Methods). When the lamellae phase is free of deformations, minimizing the free energy produces the equilibrium value where σ = d À 1 ð Þ=d, which is in agreement with complex scalar Landau models 36,37 . In our model, E is a hydrodynamic-scale field that does not involve layer spacing, so imposing an equilibrium wave number would require that covariant derivatives replace gradients [67][68][69] . At equilibrium, the free energy density is f eq = À σ In contrast to nematics, the free energy possesses two material length scales: . The coherence length ξ characterizes the defect core size and the ratio of κ = λ/ξ is a Ginzburg parameter. As in superconductors, κ<1= ffiffiffi 2 p is a type-I system, while κ>1= ffiffiffi 2 p is type-II 20 . We also take the strong anchoring limit by fixing E at solid surfaces.
Proven numerical schemes exist for minimizing the free energy of real Q tensors 70 . The numerical difficulty lies in extending the methodology to allow for complex tensor elements. We employ a gradient descent time evolution of E r, t ð Þ in 2D (see Methods 71 ). Defining the total free energy F t where μ is a mobility coefficient and Λ constrains E to be traceless and normal (see Methods). It should be stressed that E r, t ð Þ is the sole subject of all calculations-the complex amplitude ψ r, t ð Þ and layer normal N r, t ð Þ are only found ex post facto. Both |ψ| and ϕ are calculated directly from contractions of E with itself, while N is found via eigen-decomposition (see Methods). Defects are identified from the N and ϕ fields (see Methods). While our approach circumvents the ambiguity of Ψ as a double-valued function Re Ψ ½ ± iIm Ψ ½ 15,26 , the post-processing visualizations based on Re Ψ ½ introduce aberrations that E itself does not possess. It is clear that these aberrations (see Fig. 1c-e, Fig. 2 and Methods) are not visible in the |ψ| or ϕ fields, appearing only as a result of reintroducing the microscopic layer structure.

Discussion
The ability of E-theory to model simple smectics is probed by simulating a variety of confining geometries (Fig. 2). We first demonstrate that individual aspects of the order parameter can vary independently. In a long slit with N strongly anchored parallel to the walls, a temperature gradient is modeled via a linear increase of the bulk free energy (Eq. (2)) parameter A from −1 to 1 (Fig. 2a). This causes |ψ| r ð Þ to decrease, going from an ordered lamellar state on the left to the isotropic state with |ψ| = 0 on the right, in agreement with Eq. (5), and without variation of ϕ or N. After minimization of the free energy, Re Ψ ½ visualizes the lamellar structure and the layer normal reflects the direction of layering (Fig. 2a). Similarly, the phase ϕ r ð Þ can be varied within a long slit without variation of |ψ| or N by linearly increasing the phase over a narrow region of the channel walls (Fig. 2b). Again, N is strongly anchored parallel to the walls, which causes compressional distortion-the layers to dilate, as seen from Re Ψ ½ . The final pure type of deformation is distorting N r ð Þ. To create this deformation, the smectic is confined within an annulus with strong homeotropic anchoring of the layer normal (Fig. 2c). This geometry produces pure bend distortion with no compression that changes ϕ. Not all geometries allow the aspects of the order parameter to vary independently, as accomplished by Fig. 2a-c; in general, we expect some interplay between deformation modes as in Fig. 3.
There are two classes of defects in smectic systems (Fig. 1). The first relates to singularities in the layer normal N, in which N rotates by 2πm N , we show this occurring in two ways: (i) A mÑ ¼ þ1=2 disclination with a 180-degree folding of smectic layers around the defect core ( Fig. 1b; left) or a mÑ ¼ À1=2 disclination in which the layers exhibit a trifold symmetry ( Fig. 1b; center). (ii) An irregularity in the layer structure, which represents an insertion/deletion of layers at a point ( Fig. 1b; right). To identify and measure the respective charges of these, we measure winding numbers using a closed contour integral (Methods). To explore the capacity of this model to describe smectic defects, consider a circular confining domain with boundary conditions requiring a single +1/2 disclination (Fig. 1c). After minimization of the free energy, Re Ψ ½ depicts the lamellar structure around the disclination (Fig. 1f). Since N is the layer normal, it mirrors the layers (Fig. 1f, i). The lamellar structure exhibits the expected symmetries of a +1/2 disclination and deformations are primarily bend on one side of the defect and splay on the other 72 . The lamellae are highly ordered away from the defect with |ψ| ! |ψ| eq . However, |ψ| ! 0 in the defect core (Fig. 1f), verifying that E-theory permits a finite sized defect core. Crucially, E remains continuous at the center of the defect, regularizing the singularity in N. The deformations are principally curvature distortions, rather than compression, which is reflected in a constant phase everywhere in the vicinity of the disclination (Fig. 1i). We find no evidence of any artificial order parameter melting where ϕ→−ϕ, suggesting that the non-physical line tension and associated free energy penalty observed in simulations of folded layers using twodimensional scalar theories 26 is absent (Supplementary Fig. 1). The situation is analogous for a −1/2 disclination (Fig. 1d): The layers arevisualized by Re Ψ ½ (Fig. 1d), with perpendicular layer normals (Fig. 1g, j). The defect core is again seen to be locally disordered with no variation in phase, indicating negligible compression. In both ±1/2 disclinations, the free energy density is largest in the immediate vicinity of the cores (Supplementary Fig. 2). Not only is f bulk non-constant only at the core, but the deformation energy densities are strongly localized 72 . In addition to half-integer disclinations, E-theory has the capacity to simulate +1 disclinations that have been observed in transition to pairs of separated +1/2 disclinations in circular confinement with strong homeotropic anchoring ( Supplementary Fig. 3).
In addition to disclinations, which are present in both nematic and lamellar phases, edge dislocations are unique to materials with broken translational symmetry. While the phase ϕ is physically invariant to a global shift, it is set to vary linearly at the circular confining boundaries as ϕ = θ/2 for polar angle θ (Fig. 1e). This results in a dislocation: An extra layer is generated on the bottom half of Fig. 1e. While the lamellar order |ψ| still decreases in the defect core (Fig. 1h) and the order parameter variations are still localized around the core ( Supplementary Fig. 2), the phase changes by 2π around the dislocation (Fig. 1k). The occurrence of independent disclinations (Fig. 1c, d) and dislocations (Fig. 1e) highlights a strength of E-theory: since theories of ϕ alone cannot model independent disclinations and models that simulate Q near the nematic-smectic transition cannot replicate dislocations. While disclinations and dislocations are considered separately in Fig. 1, they can co-reside in a single defect.
We now consider the role of defects in lamellar states evolving to equilibrium by simulating 2D systems with a deep quench from the isotropic to lamellar state, and periodic boundary conditions (Fig. 3). At first, the system is disordered (Fig. 3a), but relaxes through defect annihilation (Fig. 3b) to form many locally ordered domains (Fig. 3c). However, even at the longest times, the system remains disordered on mesoscopic scales: It is kinetically arrested into a glassy configuration 73 with a non-zero number of defects (Fig. 3d). Within this glassy state are domains for which N is rotated by π/2 (Fig. 3c, g), which is allowed by Eq. (1) in 2D when accompanied by ϕ→ϕ±π (Fig. 3f) and correspond to the bridge-type line boundaries observed in 2D colloidal smectics 34 .
To clarify this pinning of long-lived non-equilibrium structures, we compare simulations of large and small systems. While the small system routinely relaxes to the fully ordered lamellar state (Fig. 3e; inset and Supplementary Movie 1) with lim t!1 |ψ| ! |ψ| eq , the large system never reaches the global equilibrium (Fig. 3c). Correspondingly, the free energy of the small system rapidly approaches f eq , the equilibrium defect free value; whereas, the large system is inevitably trapped away from equilibrium (Fig. 3e). Snapshots and associated videos show that both disclinations and edge dislocations are pinned 74 (Fig. 3f, g and Supplementary Movie 2). This highlights the importance of defects in lamellar ordering kinetics and the challenge posed for lamellar self-assembly [10][11][12]32,33 . In contrast to the continual relaxation Simulations initialized from an isotropic state (|ψ| ' 0 and random N and ϕ) with A = − 1 (lamellar state), C = 2, κ 2 = 0.5 and periodic boundary conditions. a-c Snapshots of Re Ψ ½ in a system size of 42ξ × 42ξ at times a t = 2μ; b 5μ; and c long-time limit kinetically arrested state (snapshot at time t = 20μ). Pink crosses (yellow trilaterals) mark + 1/2 (−1/2) disclinations. Edge dislocations with winding number ± 1 denoted by pink circles and yellow squares. d Temporal dependence of the average defect density. Shading represents the standard deviation. e Temporal dependence of the free energy density, Δf = f − f eq , for small (14ξ × 14ξ) and large (42ξ × 42ξ) systems. Shading represents the standard deviation. (inset) Steady state for 14ξ × 14ξ. f Snapshot of ϕ(r; t) corresponding to c. g The corresponding |ψ(r, t)| field. dynamics through annihilation in nematic liquid crystals 75,76 , the kinetic arrest of coarsening and long-lived domains are associated with pinned defects (Fig. 3e) 73,77 . This implies an energy barrier associated with the sliding of dislocations with respect to the lamellar structure. This indicates the possibility of non-zero Peierls-Nabarro energy barriers 26,48 , further validating of the E formulation.
The presence of inclusions embedded within the lamellar material can act to locally order layers or to induce additional defects. We evaluate boundary-induced lamellar ordering within an isotropic fluid (A > 0), due to strong anchoring to a circular inclusion (Fig. 4a). An inclusion with strong planar anchoring of N and ψ = e iπ/2 locally layers the smectic but the ordering rapidly decays (Fig. 4a). By fitting an exponential to |ψ| in a channel geometry, we extract the decay length ζ (Fig. 4b). We see that the decay length varies inversely with the Ginzburg parameter κ, indicating ζ varies linearly with lamellar coherence length ξ. While strong anchoring locally orders the isotropic phase, it induces a pair of defects in the lamellar phase (Supplementary Movie 6). The steady-state can be seen from the layer normal field (Fig. 4c) or directly from the qualitative layer visualization via Re Ψ ½ (Supplementary Fig. 4a). The topological charge of the circular inclusion is neutralized by the two −1/2 disclinations on opposite poles of the inclusion. The resemblance to a nematic system 78 follows from their shared π-rotational symmetry. Outside of the defect cores, the smectic remains well ordered and the deformation free energy contributions are localized around the inclusion ( Supplementary Fig. 5). This demonstrates the E-formalism can be employed for nontrivial geometries.
In strongly confined two-dimensional smectics composed of colloidal silica rods, it has recently been reported that half-charge ±1/2 disclinations can expand into grain-boundary lines capped by pairs of end-point ±1/4 charge defects 34 . These quarter-charge defects have been numerically reproduced using microscopic density functional theory 34 and explored using extensive Monte Carlo simulations 35,79 but are yet to be described by a mesoscopic continuum theory. To demonstrate the complex tensor E description has the capacity to simulate such defect structures, a circular inclusion is once again embedded in a simple smectic (as in Fig. 4), but initialized with a significant difference of phase ϕ at the inclusion surface compared to the bulk (Fig. 5b). This causes the −1/2 disclinations seen in Fig. 4 to each split into two −1/4 charge end-point defects capping a bridging line boundary. Such behavior is not expected nor observed in nematic liquid crystals, and does not arise in Q-tensor theories of nematics but can be reproduced by the E-tensor formalism for lamellar fluids. Both of the new quarter-charge end points are seen to possess disordered defect cores in which |ψ| ! 0 (Fig. 5a). The phase ϕ changes discontinuously across the line boundary and a π/2 misalignment of the layer normal occurs (Fig. 5b). This can also be seen in the plots of Re Ψ ½ ( Supplementary Fig. 4b). However, if the phase is integrated around the whole structure, rather than through the line boundary, the winding number of the phase of the structure is zero (Fig. 5d).
While the end-point defect cores have a free energy cost, no free energy cost is associated with the line boundary (Fig. 5c). In 2D, the eigenvalues of the traceless E differ by only a sign. This reflects the fact that a 2D smectic system can be described equally well by the layer normal N or the perpendicular layer tangent. Indeed, the layer normal and layer tangent can be swapped so long as ϕ → ϕ + π, thereby exchanging the signs of the eigenvalues. This suggests that such line boundaries are only possible in 2D, since the eigenvalue associated with the layer normal in 3D is distinct from the two degenerate eigenvalues associated with the pair of in-plane unit vectors.
Simulating the relaxation of a strongly anchored system in a annular confinements allows us to analyze the differing relaxation dynamics of both half-and quarter-charged disclinations. Systems initialized with no significant difference of phase ϕ at the boundaries compared to the bulk see formation of half-charged disclination defects (Fig. 6a, b). Annihilation dynamics are consistent across repeated initializations with noise in ψ, ϕ and N ( Fig. 6b; inset). However, systems forming pairs of quarter-charge defect complexes (Fig. 6c, d) undergo much less consistent annihilation under equally varied initial conditions ( Fig. 6d; inset), with 65% remaining after long times (t > 150μ). This is consistent with the dynamics observed in Fig. 3. By comparing the free energy density at times 5μ before and after annihilation events, the mean free energy per unit area change for a pair of half-charge defects is (5.24 ± 0.04) × 10 −3 |A| compared to (3.8 ± 0.5) × 10 −3 |A|. We note that quarter-charge defects are more costly per unit absolute topological charge than half-charge defects in these simulations, possibly accounting for their scarce creation in systems missing an enforced initial phase gradient of ϕ.
We have proposed a complex, symmetric, traceless, globally gauge invariant, normal, uniaxial tensorial order parameter E r, t ð Þ for describing simple smectic phases at mesoscopic scales. As a secondrank tensor, E encodes the apolar nature of the layer normal in an arbitrary reference frame and circumvents the ambiguity of using a complex scalar order parameter. It does so without resorting to a microscopic approach, such as density functional theory 34 , real-valued density variation 26,31 or particle-based simulations 35,79 , which also bypass such ambiguities but at a higher computational cost. While such microscopic models can simulate microscopic structure of individual layers, a numerically amenable framework for simulating the mesoscopic variations is advantageous for modeling configurations with layer normal shown in red. b Inverse of the exponential decay length ζ as a function of Ginzburg parameter κ for the isotropic phase confined between planar walls. c Plot of |ψ| in the vicinity of an inclusion (same parameters as a except R = 4ξ and A = −1 (lamellar state)). The system is initialized near isotropic (|ψ| ∼ 0:2). and dynamics at large length and time scales. By conjoining local layer orientation and the extent of ordering into a single mathematical object, the E tensor is a mesoscopic description that can reproduce both disclination and dislocation defects with finite defect cores.
Though individual singularities can be elegantly analytically handled through local branch cuts, the tensor order parameter description globally eliminates this ambiguity in a numerically pragmatic manner. Akin to the nematic Q tensor, this has the numerical advantage of avoiding point singularities. In the presence of disclinations, the layer normal N is a singular vector, while the phase ϕ is singular about dislocations-the tensor E regularizes the singularities and remains continuous at the center of defects. In this paper, we have restricted our consideration to a "simple" smectic (purely lamellar) idealization that does not require a well-defined nematic director as a prerequisite of spontaneous lamellar symmetry breaking. Liquid crystalline phases such as smectic-A or -C are composed of anisotropic molecules that also exhibit nematic ordering, and simulating these will require coupling of E to Q, which is conceptually straightforward 36,80 .
Additionally, the E-formalism should be extended to consider more complex systems. Smectic textures in three dimensions can be complex 81 , appealing 82 and difficult to model 31 . The tensor theory could be extended to simulate patterned defect arrays through programmable photoalignment 83 , electrically reversible templating 84 and micropatterned substrates 18 . Much like the introduction of Q helped expand the possibilities for numerical modeling of nematics, we expect this framework to be advantageous for simulating colloidal smectics 82,85-90 , smectic-isotropic interfaces 91,92 , smectic-smectic emulsions 93 , smectics in contact with active material 94 and swimming bacteria in smectics 95 .

Numerical methods
We employ a gradient descent evolution of E r; t ð Þ. This is to say that E r; t ð Þobeys a time-dependent Ginzburg-Landau model, which follows the steepest decrease in global free energy under the constraints that it remain traceless and uniaxial. This is described by Eq. (6) in the main text and Eq. (16) in Methods.
The total free energy density is the sum of Eq. (2)-(4), At each timestep, we use the instantaneous free energy density to calculate the right hand side of Eq. (6). The Lagrange multipliers can be found directly but the functional derivative term involves spatial derivatives. The system is discretised in space on a square grid of step size Δx. We employ a two-step Adams-Bashforth method to calculate E at the next timestep, using a discrete time step of Δt. This is iterated for T time steps. For Fig. 2, 3, Δx = 0.7ξ, Δt = 0.001μ. For Fig. 2, T = 15000, which amounts to a total simulation time of 15μ; for Fig. 1, T = 20,000 We employ specific boundary and initial conditions for each system considered, as discussed in the main text. To initialize the system, we specify E at each point in space, which we do by setting a local ϕ, |ψ| and N. These are then converted to E through the definition Eq. (1).
After initialization, we directly simulate E, without any reference to ϕ, |ψ| and N. Only after the simulations have completed are ϕ, |ψ| and N found from the resulting E. The complex scalar order parameter ψ corresponds to one eigenvalue and N is the associated eigenvalue. Due to the tracelessness of E, the pair of eigenvalues contain the same information and differ only by their sign in 2D. In 3D, the pair of inplane eigenvectors (i.e., the pair that are orthogonal to the layer normal N) must be interchangeable and so their associated eigenvalues are degenerate, while the third is constrained to be proportional by the tracelessness condition.
As complex eigen-decomposition is numerically imprecise, we find the following method more reliable. The scalar order parameter |ψ| is found directly by contracting the complex tensor with its complex conjugate, E : E * = σ|ψ| 2 . Explicitly, |ψ| = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E ij E * ji =σ q . Likewise, the phase ϕ is found by contracting the complex tensor with itself E : E = σ|ψ| 2 e 2iϕ , o n ce |ψ| is known. Explicitly, ϕ = ðÀi=2Þ arg½E ij E ji =ðσ|ψ| 2 Þ, which returns ϕ ∈ [−π, π). Finally, to determine the layer normal direction N, we apply eigen-decomposition to the real tensor e −ϕ E, recovering a pair of real eigenvectors. In the standard case, the eigenvector corresponding to the positive eigenvalue is selected as N. However in systems with boundary conditions enforcing a phase discontinuity, such as in Fig. 1c, we must enforce the choice of eigenvalue along each boundary and promulgate that numerical decision though the system. This decision prevents an artifical discontinuity in the layer normal along an arbitrary line due to a phase shift of π swapping the signs of the eigenvalues, and hence the choice of corresponding eigenvector. This argument is avoided in 3D as the eigenvalues are ψ, −ψ/2, −ψ/2; so the eigenvalue with the largest modulus should always be chosen.

Lagrange multipliers
We enforced the fact that E be traceless, uniaxial and a normal operator, which commutes with its own adjoint, in the numerics. This is essential for ensuring E maintains real eigenvectors and allows us to interpret the E after the simulations. These two conditions can be written as where [A, B] = AB − BA denotes the commutator. Using the Cayley-Hamilton Theorem and noting g 2 can be rewritten Decomposing g 1 into real and imaginary parts, g 1a = Re g 1 Â Ã and g 1b = Im g 1 Â Ã , we introduce three real Lagrange multipliers, λ 1a , λ 1b and λ 2 , defined as Where c 1 and c 2 are and we have written E = X + iY for real tensors X and Y. We then introduce a dynamics that minimizes our free energy above with respect to these constraints With respect to the main text, Comparison to existing models Substituting the eigenvalue ψ and vector N into Eq. (1) produces explicit forms for the free energy densities in terms of |ψ|, ∇ϕ and N. This allows us to compare to existing models. The bulk term (Eq. (2)) becomes where σ = d À 1 ð Þ=d, which demonstrates the consistency between this complex tensor theory approach and scalar-based bulk free energies 67 . Other possible real terms of the form ∼ trðE + E * Þ β are zero, while terms of the form ∼ ½E + E * : ½E + E * depend directly on ϕ, which is nonphysical. Therefore, such terms cannot be included in Eq. (2). The compression term (Eq. (3)) becomes where we have make use of the identity N i ∇ j N i = 0, which is due to the fact that N is a unit vector. Differentiating this again, we see that , which we make use of below. Finally, the curvature term (Eq. (4)) is the most complicated, becoming where we have stated the contraction on the gradients in the last term using Einstein notation for clarity but have used vector notation elsewhere. It is worth re-emphasizing that none of the contributions to the free energy density depend directly on the phase ϕ. In these forms, the E-tensor theory can be compared to existing models. In the ground state equilibrium of flat, equally spaced, layers, the deformation free energies (Eq. (19) and Eq. (20)) are both zero. Only the bulk free energy f bulk (Eq. (18)) is non-zero and its form is consistent with scalar-based bulk free energies 36,37,67 .
In the limit of fixed ψ but variable N, only incompressible distortions are allowed. Taking ψ to be constant, the deformation free energies (Eq. (19) and Eq. (20)) become If we further only keep the lowest order term in the deformation free energy ∼ ∇N ð Þ : ∇N ð Þ, then this is precisely the nematic deformation free energy in the one-constant approximation 70 . We note that this specifically does not lead to the Oseen constraint that twist be prohibited, a prevalent simplifying assumption in models of smectic materials, because we have explicitly made a one-constant elastic approximation for simplicity. In nematic models near the vicinity of the nematic-smectic A phase transition, the ratios of elastic constants significantly differ from unity 96 . More intricate E-tensor theories that allow for differing elasticities will be able to make twist and bend of the layer normal (splaying of the layers themselves) come at a increased free energy cost.

Biaxiality
In the present study, we have exclusively focused on 2D systems for which uniaxiality is maintained due to the tracelessness constraint (requiring eigenvalues to be equal and opposite). However, a degree of biaxiality could be conceived along a secondary direction in 3D. In three dimensional nematic liquid crystals, allowing biaxiality increases the number of degrees of freedom in the order parameter Q from three in the uniaxial case to five. The two ancillary degrees of freedom represent the degree of biaxial alignment and the biaxial direction, constrained to be both a unit vector and orthogonal to the director. However, bulk biaxiality at equilibrium requires higher order terms in the expansion of the bulk free energy 64 . On the other hand, the higher order terms are not needed for biaxiality to naturally emerge in nematic defect core regions 108 . This is because the three eigenvalues can be distinct in regions where the ordering goes to zero 109 and in strong confinements 110 . Similarly, biaxiality in E increases the degrees of freedom from four in the uniaxial case to seven. The three extra degrees of freedom represent another degree of ordering, biaxial layer displacement and the secondary direction, constrained to be a unit vector orthogonal to N. The biaxial case will not necessarily have a uniform complex phase across different components of the tensor and biaxiality at equilibrium would require higher order terms in the bulk free energy expansion (Eq. (2)).

Defect identification
To identify defects, we measure winding numbers at each point over the smallest possible closed loop Γ, as a measure of topological charge. For disclinations, we measure the change in the layer normal azimuthal angle, dα = N ⋅ dl. Values of m N = ± 1/2 correspond to half-charge topological defects. For dislocations, α = ϕ measures the displacement of layers and gives the strength of the Burgers vector.

Layer visualization
To visualize the layers, we use Re Ψ ½ = |ψ| cosðq 0 N Á r + ϕÞ: ð26Þ Note however that this visualization technique does not respect the smectic symmetry; the focus is on providing an intuitive visual description of the layer structure. To calculate Re Ψ ½ , we use a Voronoi transformation of the plane into regions by defects and use the location of each defect as the origin for the dot product in Eq. (26). This numerical scheme gives good visualization of isolated defects in terms of layer structure (see Fig. 1c-e) but the transition from smooth |ψ| and ϕ fields can introduce aberrations in defect-crowed regimes [see Fig. 3c]. This is due to the ambiguity of the layer displacement that is not possessed when working in terms of E alone. Visualizations based on Re Ψ ½ must necessarily re-introduce the shortcomings that E-theory bypasses and so we generally avoid visualizations of Re Ψ ½ , except as a qualitative guide. Despite these disadvantages, we at times find it convenient to have an explicit, if rough, visualization of the layer configuration. The advantage of working with these fields can be seen in Fig. 3f, g and Supplementary Movies 2, 3. The resulting layer visualization is in Fig. 3c and Supplementary Movie 4.

Data availability
All data generated in this study are included in this published article (and its supplementary information files).

Code availability
Codes are available upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.